Mechanism of murine epidermal maintenance: Cell division and the Voter Model 
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The dynamics of a genetically-labelled cell population may be used to infer the laws of cell 
division in mammalian tissue. Recently, we showed that in mouse tail-skin, where proliferating 
cells are confined to a two-dimensional layer, cells proliferate and differentiate according to a simple 
stochastic model of cell division involving just one type of proliferating cell that may divide both 
symmetrically and asymmetrically. Curiously, these simple rules provide excellent predictions of the 
cell population dynamics without having to address their spatial distribution. Yet, if the spatial 
behaviour of cells is addressed by allowing cells to diffuse at random, one deduces that density 
fluctuations destroy tissue confluence, implying some hidden degree of spatial regulation in the 
physical system. To infer the mechanism of spatial regulation, we consider a two-dimensional model 
of cell fate that preserves the overall population dynamics. By identifying the resulting behaviour 
with a three-species variation of the "Voter" model, we predict that proliferating cells in the basal 
layer should cluster. Analysis of empirical correlations of cells stained for proliferation activity 
confirms that the expected clustering behaviour is indeed seen in nature. As well as explaining how 
cells maintain a uniform two-dimensional density, these findings present an interesting experimental 
example of voter-model statistics in biology. 

PACS numbers: 



I. INTRODUCTION 

A major challenge in biology is to determine how pro- 
liferating cells behave in developing and adult tissues. To 
gain insight into the processes of cancer onset, aging and 
wound healing, biologists have long recognised that the 
spatial organisation of cells in tissue provides indirect ac- 
cess to the underlying cell behaviour. Some tissues, such 
as the auditory hair cells of the inner ear are arranged 
into repeating units containing groups of specialised cells 
essential for the function of the tissue pQ. In contrast, 
in other tissues cells do not organize into coherent struc- 
tures that reflect their cooperative function. Indeed, the 
arrangement of some cell types appears random [2]. In- 
ferring the rules of cell behaviour in these apparently 
unstructured tissues appears challenging. One may ask, 
therefore, how cell behaviour in such tissues is regulated 
in the absence of well defined spatial roles. 

In this context, it is interesting to consider the case of 
mouse tail-skin epidermis, for which the laws governing 
cell behaviour have recently been resolved. Mammalian 
epidermis comprises hair follicles interspersed with in- 
tcrfollicular epidermis (IFE), which consists of sheets of 
specialised cells known as keratinocytes [3J, see fig.[lja). 
Cells are shed continually from the epidermal surface, 
and are replaced by proliferating cells in the basal layer, 
whose progeny may cease proliferating and then migrate 
through the suprabasal layers before reaching the epi- 
dermal surface. Until recently, it was widely assumed 
that adult tissue is maintained by two different prolifer- 
ating cell populations. These comprised long-lived, self- 
renewing stem cells, that have the potential to undergo 
an unlimited number of cell divisions and which main- 
tain a second population of transit-amplifying cells (TA), 
whose proliferative potential is limited [4 . After several 



rounds of division, it was conjectured that TA cells dif- 
ferentiate, exit the cell cycle and move out of the basal 
layer. 

In a recent study by Clayton et al, inducible genetic 
labelling was used to study the mechanism of epidermal 
maintenance by tracking the fate of a representative sam- 
ple of cells and their progeny (clones) in normal murine 
tail epidermis [5]. By analysing the size distribution of 
such clones over a period of one year, it was possible to 
infer that the epidermis is maintained by just one type 
of progenitor cell. According to the revised model of 
epidermal maintenance, progenitor cells capable of both 
symmetric and asymmetric division give rise to a popu- 
lation of non-proliferating cells, which then transfer from 
the basal layer to the suprabasal layers. More precisely, 
labelling proliferating cells as type A and differentiated 
basal layer cells as type B, the basal layer cell population 
is governed by the stochastic non-equilibrium (Marko- 
vian) process, 

!A + A Prob. r 
A + B Prob. 1 - 2r 
B + B Prob. r M 

B , 

involving three adjustable parameters: the overall cell 
division rate, A, the proportion of cell divisions that are 
symmetric, 2 x r, and the rate of transfer, T, of non- 
proliferating cells from the basal to the suprabasal lay- 
ers. As well as overturning the accepted paradigm of 
epidermal homeostasis being achieved by discrete popu- 
lations of stem and TA cells, the model provides a de- 
gree of quantitative predictive rigour that is unusual in 
the field of cell tissue biology; for example by quanti- 
fying the division and migration rates (A = 1.1 /week, 
r = 0.31/week) and the branching ratio (r = 0.08), as 
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FIG. 1: (Color online) (a) Schematic cross-section of murine 
interfollicular epidermis (IFE). Proliferating cells (grey) are 
confined to the basal layer (labelled i); differentiated cells mi- 
grate through the superbasal layers (ii), where they flatten 
into cornified cells, losing their nuclei and assembling a corni- 
fied envelope (green online) (iii), eventually becoming shed at 
the surface, (b) Two examples of typical clones acquired at 
a late time point, viewed from the basal layer surface. Cell 
nuclei are labelled blue; the hereditary clone marker (EYFP) 
appears yellow. Scale bar: 20/im. 



shown in Refs. 015]. As a result of its quantitative na- 
ture, the model establishes a platform for investigating 
the role of different cellular constituents (such as gene 
products) in regulating cell fate, for example by study- 
ing how different genes influence the model parameters. 

It is an interesting fact that process is capable 
of fitting the wide range of clone fate data within a 
"zero-dimensional" framework, i.e. without having to 
address the spatial orientation of cells within the basal 
layer. Yet, the observed uniformity of cell density im- 
plies a degree of regulation beyond that which can be 
addressed in the zero dimensional framework. In par- 
ticular, when augmented by spatial diffusion, the pro- 
posed model leads to "cluster" formation in the two- 
dimensional system whereupon local cell densities are 



predicted to diverge logarithmically [7J [H] . In biological 
terms, such behaviour would correspond to a severe dis- 
ruption of the epidermis, with much of the tissue dying 
away, leaving only a few isolated and very thick clusters 
of epidermis. Significantly, the observation that labelled 
families of cells remain largely cohesive (see for exam- 
ple Fig. T|d), reveals that cell mobility must be small, so 
that such divergences would be significant within a mam- 
malian lifetime. These divergences can not be regulated 
through a local density-dependent mobility. 

Thus, the success of the zero-dimensional fit, despite 
the predicted divergence in two dimensions, leaves us 
with an interesting challenge that is the focus of the 
present study: Can we uncover, from the spatial dis- 
tribution of basal layer cells, the mechanism by which 
cells regulate a uniform cell density without compromis- 
ing the integrity of the zero-dimensional fit, as embodied 
in process ([T])? 

In order to identify the underlying rules of cell division 
and differentiation, we shall draw upon the results of two 
types of experiment. First, we shall revisit the clone fate 
data used by Clayton et al, and examine the previously- 
discarded spatial distribution of labelled basal layer cells 
for signatures of underlying regulation. Second, we shall 
consider the statistics of the entire population of basal 
layer cells. In particular, by immunostaining basal layer 
cells for markers of cell proliferation, it is possible to 
analyse the spatial distribution of all progenitor cells. 

Thus, the aim of this paper is to elucidate how the 
experimental observations constrain any proposed the- 
ory of spatial behaviour in the basal layer. In summary, 
we shall show that the dynamics predicted by process 
are indeed consistent with the constraint on uniform 
cell density, provided that cell division occurs only upon 
the migration of a nearby type B (i.e. differentiated) 
cell into the basal layer. Moreover, we confirm that the 
clone fate data is consistent with a restricted degree of 
cell mobility, whereby cell motion is not diffusive and 
random. Instead, differentiated cells only migrate lat- 
erally as a response to fluctuations in the local density. 
Finally, to test the validity of the proposed spatial model, 
we use it to predict that, while maintaining a uniform 
total areal cell density, the population of progenitor cells 
should cluster over time. By considering the radial corre- 
lation function for the spatial distribution of progenitor 
cells, we find that this prediction is in good qualitative 
agreement with experiment. Quantitatively, the compar- 
ison reveals that the experimental degree of progenitor 
cell clustering is slightly higher than that expected for 
the parameter value of r = 0.08 determined previously 
through clonal analysis. Although several technical diffi- 
culties may challenge the reliability of these quantitative 
results, we speculate that such excess clustering may be 
a signature of spatial regulation of cell fate during asym- 
metric division. These results also shed light on previous 
observations of clustering of cells undergoing mitosis in 
the epidermis [HI [TO] ■ Such observations have been in- 
terpreted in the biological community to be a signature 
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FIG. 2: (Color online) Schematic motivating the proposed 
lattice model defined in Eq. Q. Top: A cartoon of cells 
within the basal layer, showing the exit of a cell (light grey) 
through upward migration into the suprabasal layers, concur- 
rent with the division of a nearby progenitor cell (red online) . 
To ensure continuity of the basal layer, it is postulated that 
cells rearrange to maintain a uniform cell density (wide ar- 
row). Bottom: In a simple lattice model that captures the 
essence of the steady-state dynamics, the exit of a cell from 
the basal layer gives rise to a vacant lattice site (light grey), 
which rapidly diffuses by exchanging position with adjacent 
cells (white). Upon coming into contact with a proliferating 
cell, the latter may divide and replace the vacancy with a 
daughter cell (striped). 



of regulation that leads to coordinated cell division. By 
contrast, this work shows that the tendency of prolif- 
erating cells (and therefore mitoses) to cluster is in fact 
consistent with cells dividing independently and stochas- 
tically — indeed it is the hallmark of the proposed spatial 
process. 

This paper is organised as follows. In section |II A| 
we develop a phenomenological model of cell behaviour 
that incorporates the experimental constraint on uni- 
form cell density. We identify the proposed model 
as a variation upon the so-called monomer-monomer 
model of surface catalysis. We then analyse this sim- 



pler model in section II B including an exact solution for 
the two-point density correlations in the closely-related 
monomer-monomer model. In section IIIII we test the 
uniform-density model against a range of experimental 
data, first through a qualitative comparison of the model 
with the empirical clone shape data (section |lll A ), and 



then by analysing the spatial correlations observed for 
proliferating cells (section IIIB I. We conclude with a dis- 



cussion of cell clustering in section IV 



II. A SPATIAL MODEL OF CELL KINETICS 

To understand how the experimental observations con- 
strain any proposed theory of spatial behaviour in the 



basal layer, we shall first address the constraint imposed 
by the observed uniform cell density. To this end, we 
postulate that proliferating cells only divide upon the 
migration of a nearby differentiated cell into the supra- 
basal layers (see Fig. [2]). This requirement is purely phe- 
nomenological, as it ensures a uniform density without 
specifying the mechanism by which it is implemented. 
Indeed, a range of regulatory pathways can be seen to 
give rise to the same phenomenology, for example by cou- 
pling cell division processes to the local stress [TTJ [12] or 
using short-ranged morphogen gradients as well as feed- 
back by cell-cell communication. Once this initial con- 
straint is accounted for in this section, we shall draw 
upon the further observations of clone and basal layer 
morphology in order to identify additional rules govern- 
ing cell behaviour. 

In general, it is a challenge to couple the division 
and migration of basal layer cells while still allowing for 
some degree of cell compressibility, whereby a dividing 
cell may compensate for the exit of a non-adjacent cell 
through lateral motion (see Fig. [2] top). Several ap- 
proaches have been used in the past to overcome this 
problem. In the context of tissue development, one may 
treat the cell tissue as an elastic medium, whereby the 
local cell density is coupled to the cell division process 
through the stress in the surrounding tissue [IT]. Such 
an approach is capable of accounting for a range of re- 
alistic properties of two-dimensional cell tissue growth, 
such as the distribution of cell sizes, as well as cell com- 
pressibility. Yet, for the simple problem of steady-state 
tissue maintenance, involving no net growth, it is un- 
likely that the complexity of the elastic tissue model is 
required to explain much of the experimental data. A 
second approach treats the two-dimensional basal layer 
as a 'foam' of cells (a Voronoi tcssclation) , for which the 
steady-state condition may be used to relate between 
the local cell topology (or more precisely, the number 
of nearest-neighbours in the basal layer) and the likeli- 
hood of division or migration. Applied to the problem of 
epidermal maintenance, this approach successfully pre- 
dicts the steady-state topology of epidermal basal-layer 
cells [13] , and it identifies that cells with a larger number 
of neighbours are more likely undergo division, whereas 
cells with a small number of neighbours are more likely 
to migrate in to the supra-basal layers. However, it is a 
challenge to extend this model to allow for two distinct 
cell populations that are exclusively committed to cither 
division or migration, such as described by process ([I]). 
Yet a third approach draws upon simulations in which 
cells are modeled as quasi-spherical particles that de- 
form during cell division [14] . However, as this approach 
draws upon a wide range of (uncontrolled) parameters 
to describe the cell-cell interactions, it is more complex 
than required for this case. 

Therefore, in the following we shall suffice with a sim- 
pler description of the basal layer, by drawing upon the 
non-equilibrium lattice models discussed in the recent 
literature. In particular, we shall model the basal layer 
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as a lattice in which each site is occupied by one cell. 
Cell compressibility is then modeled by a population of 
lattice vacancies, which are created upon cell migration 
(B — > 0) and then diffuse rapidly as compared to the 
cell kinetic rates (A, T), before annihilating upon the 
division of an adjacent cell (see Fig. |2|. Since this rule- 
based model is not primarily based on a direct physical 
representation of individual cells, it may overlook cer- 
tain physical effects. For example, the migration of type 
B cells out of the basal layer may be facilitated by me- 
chanical forces exerted by neighbouring cells — a situa- 
tion that is hard to represent with a cellular automaton. 
However, recalling that the stochastic rules embodied by 
process Q have been experimentally verified, it is rea- 
sonable to start by considering a similar stochastic pro- 
cess in two dimensions. Later, we shall further justify 
the use of the lattice model by showing that the lattice 
geometry does not effect the qualitative behaviour of the 
system (see section IlfTe)). 

Although a significant advantage of the lattice descrip- 
tion of the basal layer is the ease by which it may be sim- 
ulated, a range of analytic results are also made accessi- 
ble by showing that in certain limits, the model reduces 
to a simple two-component model that belongs to the 
generalised voter model universality class [15j [THJ ITF] . 
The voter model universality class describes lattice pro- 
cesses that undergo phase separation in two dimensions 
in the absence of surface tension. In the context of basal 
layer kinetics, this "phase separation" corresponds to 
the clustering of proliferating cells. Curiously, in the 
special case r = 1/4, the basal layer lattice model re- 
duces to the reaction-limited monomer-monomer model 
of surface catalytic reactions, in which the classical voter 
model dynamics are augmented by infinite-temperature 
Kawasaki exchange dynamics [IB]. As well as drawing 
upon a range of existing results afforded by these models, 
we shall derive here an exact solution for the two-point 
spatial correlation function of the monomer-monomer 
model. These results reveal the continuous transition 
between voter-like behaviour and diffusive behaviour as 
the relative rates of symmetric vs. asymmetric cell di- 
vision (or reactant deposition in surface catalysis) are 
adjusted. 



A. The lattice model 

As mentioned above, we are interested in constructing 
a spatial model that recovers the behaviour of process 
, and maintains a uniform cell density. To account for 
the steric repulsion of basal layer cells, we will charac- 
terise the basal layer as a two-dimensional lattice, where 
each site is host to one of the two cell types, or it remains 
vacant. When the vacancy fraction is very low, then such 
a lattice description presents a reasonable approximation 
of the observed near-uniform arrangement of basal layer 
cells during normal adult skin maintenance. Then, to 
regulate the cell density, progenitor cells (A) are allowed 



to divide only when neighbouring a site vacancy, while 
the migration of post-mitotic (B) cells from the basal 
layer leads to the creation of vacancies which are free 
to diffuse in the basal layer through the displacement of 
neighbouring cells (see Fig. j2|. In summary, denoting a 
site vacancy with the symbol 0, the lattice model may 
be written in terms of the non-equilibrium process: 



A A Prob. r 



A 
D 

X 




(2) 



where the site hopping rate a reflects the capacity of 
vacancies to diffuse within the basal layer, and X is used 
to denote either a type A or a type B cell. To gain 
some initial insight into the dynamics, and to identify 
constraints on the parameter space, one may consider 
the steady-state mean-field cell densities associated with 
this process. It is straightforward to show that the mean- 
field equation for the the vacancy fraction n is given by 



d t n = o-V 2 



r(i 



)) - A'pn , (3) 



where p is the (constant) A-cell fraction. From here one 
may sec that the near-uniform cell density in the basal 
layer (n <C 1 and uniformly distributed) constrains us 
to the region of parameter space A' ^> T, such that any 
vacancy created through the migration of a type B cell 
out of the basal layer is rapidly removed upon coming 
into contact with a proliferating cell. In this limit, the 
numerical value of the parameter A' becomes irrelevant, 
as may be seen by calculating the effective local division 
rate A = A'n , which takes the value A = T(l — p)/p, 
independent of A'. This relationship between the rate of 
cell division and migration is identical to that obtained 
in process (JTJ [B]. 

Although the uniform vacancy density is a stable fixed 
point of the mean-field dynamics (Eq. [3]) , one may worry 
whether fluctuations about the mean-field solution are 
capable of compromising the uniform density of the basal 
layer cell lattice in actual manifestations of process Q . 
To eliminate this concern, one requires rapid dissipation 
of density fluctuations independently of the cell kinetics, 
from which we infer a that the lattice vacancy population 
must diffuse rapidly compared to the time-scale of cell 
division and upward migration, viz. a T, as stated 
earlier. Biologically, this condition corresponds to the 
assumption that cells are largely incompressible, so that 
local density fluctuations lead to the rearrangement of 
cells on a time-scale that is significantly faster than that 
of cell division (cf. Ref. [H]). In section III A we shall 



show that the empirical clone fate data further constrain 
the parameter space to the region A' 3> c, whereby the 
exit of a type B cell is compensated for by a nearby cell 
division. 
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With these definitions, in the parameter space A' ^> 
a r one may see that the spatial model introduces no 
new relevant parameters compared to process pi. That 
is, the model behaviour depends only on the (known) 
zero-dimensional parameters (r, A, T), with the contri- 
bution of the new parameters (A', a) entering through 
the dimensionless ratios T/X', T/a and u/X' , all of which 
vanish. But, as a precondition on its validity, does the 
model also reproduce the observed 'zero-dimensional' 
basal layer clone size distributions? This is by no means 
obvious, given the critical (and therefore delicate) na- 
ture of process ([TJ [B]. For example, a high density of 
progenitor cells in the lattice model may lead to vacancy 
depletion and jamming, an effect that has no analog in 
the zero-dimensional system. Therefore, in section |III| 
we shall show, using Monte-Carlo simulations, that the 
proposed lattice model indeed succeeds in reproducing 
the empirical clonal statistics from Ref. [5]. With this 
basic confidence in the validity of the model, we now 
proceed to analyse its behaviour in more detail. 

In general, the cell kinetics in process Q describe 
a hard-core non-equilibrium system involving three cell 
species. Recent progress in the field of non-equilibrium 
statistical mechanics has resulted in several possible for- 
malisms with which to study such systems [T9j |2Ql El] . 
However, for the case at hand, these approaches are 
unnecessarily complex. Rather than analysing the cur- 
rent microscopic model, it is convenient to recast the 
cell kinetics into a simpler form that describes the same 
phenomenology In particular, for the parameter space 
of interest (A' 3> o T), it is sufficient to consider 
a lattice fully occupied by A and B cells, without ad- 
dressing the population of lattice vacancies. To see 
this, one may see from Eq. (|3| that the vacancy dy- 
namics occur on a fast time scale compared to that of 
cell migration (1/r), and division (1/A). Therefore, re- 
ferring back to the lattice process we may (heuris- 
tically) eliminate the vacancy population by replacing 
the cell division process (A —> X Y) with a direct 
cell-cell 'reaction' process (A B — > X Y), and by re- 
placing the A-cell division rate A' with the effective rate 
A(x) = A'n (x) = r[l — p(x)]/p(x). Here, the local 
A-cell fraction p(x) refers to the A-cell number density 
coarse-grained over the nearby lattice neighbourhood. 
For example, denoting the number of type A cells on 
a lattice site as n x , we have p(x) = n x /u>(|x — x'|) 
where w(x) is some suitably chosen normalised envelope 
function. Within this framework, one may then replace 
process ^ with the more simplified form 



A A Prob. r 

. A(x) I A B Prob. | -r 

B A Prob. l-r 

B B Prob. r 



(4) 



The degree to which this heuristic simplification is indeed 
justified will be discussed at the end of section |IIB| to- 
gether with a quantitative comparison of the behaviour 
of the exact and simplified models. It is already clear 



that process Q cannot describe the explicit upward mi- 
gration of post-mitotic cells from the basal layer. How- 
ever, in the physically relevant limit A' 3> c, eliminating 
the vacancy population has no qualitative effect on the 
statistics of the progenitor cell compartment, and there- 
fore processes ^ and Q are expected to result in the 
same basal layer phenomenology. 

Interestingly, process Q is closely related to the model 
of monomer-monomer surface catalysis [22 , 23J . In par- 
ticular, when the coarse-grained distribution of type A 
cells is effectively uniform, such that A(x) ~ const., then 
one may identify the symmetric branches of process Q 
with the classical zero-temperature voter model, while 
the asymmetric division channel AB^BA, describes 
the Kawaski dynamics of an infinite-temperature Ising 
spin model. Later it will become clear that for the em- 
pirical value of r — 0.08, the significant contribution of 
the latter will justify the approximation of near-constant 
A. 

This analogy provides access to several known results. 
In Ref. |18j . Krapivsky showed that starting from ran- 
dom initial conditions, the classical voter model (i.e. 
with r = 1/2) will lead to a lattice of N sites becoming 
completely saturated with either type A or type B cells 
after a time t ~ NlnN. Moreover, Frachebourg and 
Krapivski gave an exact solution for the two-point spatial 
correlations in this case [22], from which they inferred 
that, in the time leading up to saturation (t <C NlnN), 
the different cell types separate into domains of ever- 
increasing size, with a typical lengthscale growing as 
L ~ lni, and with the density of interfaces cab between 
type A and type B cells dropping as cab ~ 1/ In i. As the 
system approaches saturation (t ~ N In TV) , one of the 
cell types comes to dominate. The classical voter model 
is an example of domain growth in the absence of surface 
tension [23]. Therefore, the boundaries between domains 
rich in A and B cells are completely unstable, leading to 
strikingly different and irregular domain morphologies, 
compared to the smooth phase-separated shapes result- 
ing from surface-tension mediated domain growth. 

Qualitatively, the results found for the classical voter 
model (r = 1/2) allow us to make several interesting 
predictions relating to the spatial distribution of A and 
B cells. In particular, some degree of clustering of pro- 
liferating cells is to be expected in adult mice, result- 
ing from the growth of domains rich in progenitor cells. 
Moreover, the ongoing growth of the domain size L sug- 
gests that larger clusters are expected in old vs. young 
epidermis. Yet, to make full contact between process ^ 
and the empirical data, it becomes necessary to calculate 
the model properties whilst allowing for the relatively 
low value of r = 0.08 found in the experimental system. 
Therefore, in the following, we shall extend the analysis 
of Frachebourg and Krapivsky to obtain results valid for 
arbitrary r. Indeed, with r = 1 /4, the following analy- 
sis results in the exact solution to the reaction-limited 
monomer-monomer surface catalysis model. 
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B. Exact solution for two-point correlations 

In the following, we will follow the same approach as 
taken in Refs. [HI E2] for the r = 1/2 case, but we 
generalise to allow for arbitrary r and different lattice 
geometries. For completeness, we include here aspects 
of the solution that were also described in some detail in 
Ref. [IH] , such as the Master equation and the dynamical 



equations required to define the problem. We start by 
identifying type A cells with state 1 and type B cells 
with state 0, so that a lattice with site index i may be 
described in terms of the Ising variables $ = {rti}, rii € 
{0, 1}. Referring to Ref. [TS], the Master equation for the 
probability distribution P(3>, t) for the system to occupy 
state <I> at time t is given by, 



ul e) (F^)P(F^,t) + U± e >(F l+e $)P(F l+e <f>,t) 



<\ - r)ul e) (F i F i+e 'i>)P(F i F i+e ^,t) - (~ + r)uf\<Z>)P($,t) 



(5) 



Here, {e} represent the nearest-neighbour lattice vectors 

(|e| = 1), and J// e) ($) <E {0,1} indicates whether the 
cells at sites i and i + e are a 'reactive' pair, viz. 

f// e) ($) = m + m+e - 2nim+e . 

The "spin-flip" operator is defined by F^ — 
{rij for all j ^ i; 1 — n^, so that F^ and F i+e <& corre- 
spond to the symmetric division channels and FjFi +e <& 
corresponds to an asymmetric division in which the lo- 
cation of the type A and B cells is reversed (viz. AB^ 
B A). 

From here, recalling that (rairij) = nirijP(&, t), it 
is simple to show that the two-site correlation function 
evolves according to the discretizcd diffusion equation 
for non- neighbouring sites, 

j t (n i n i ) = ^(A i + A i )(n i n i ), (6) 

where Ai is the discrete Laplacian operator, defined by 
AiUi = X) e ( n i+e — ft*0- However, the diffusion equation 
is modified for nearest-neighbour correlations, giving 

^(?M«i+e) = ^(A ; + A i+ e)(nirii +e ) 

+ (o _ r)\[2(riin i+e ) - (TH+e) 

<m>], (7) 

and the on-site moment is trivially (njrii) = (ni). 

Making the simplifying assumption that the initial dis- 
tribution P(<1>,0) is translationally invariant, then the 
fraction of type A cells is given by (ni) = p = const., 
and (ninj) depends on (i — j) at all times. Therefore, 
introducing the correlation function C\ — (njTjj+i), we 
can rewrite Eqs. ^ and Q as follows: 

^Ci - AA ; Ci - Ml - 2r)A(p - C e ) , (8) 



for |i| > 1, subject to the constraint Co = p = const. 
That is, the correlation function evolves according to 
a discrete diffusion equation with sink terms at the 
nearest-neighbour sites and with a fixed boundary con- 
dition at the origin. The linear nature of the prob- 
lem allows one to seek a solution in terms of the rel- 
evant Green's function, e.g. Ar 1 = \(t) = 
e~ t Ii tB (2Xt)Ii y (2Xt) for a square lattice i = (i x ,i y ). 
For uncorrelated initial conditions, viz. Ci(t = 0) = 
p 2 + 5i a(p — p 2 ) , one may write down a general solution 
in the form 

Ci(i) = p 2 + { P -p 2 )G i {t) 

+ E / *-Jj(r)G,_j(t-r), (9) 
. Jo 

where Ji(i) is the source distribution required to both 
maintain Co = const, and also to incorporate the sink 
terms from Eq. ds]) . Thus, J t (t) = for |i| > 1, J e = 
— (1 — 2r)X(p — Ce), and J — zX(p — C e ), where z is the 
number of nearest neighbours, and the initial conditions 
imply that C e is the same for all nearest neighbours. For 
now we will explicity consider the square lattice [z = 
4), however the effect of changing lattice geometry will 
be discussed near the end of this section. Making use 
of Eq. (§, one may write down a set of self-consistent 
equations for the source terms, viz. 

J (t) = AX^p-p 2 -J^ [j (t-r)G (04) (r) (10) 
+ J e (t - t) (G (0i0 )(t) + G (0 , 2) (t) + 2G (M) (r)) ] j 

Mt) - - 1 ^j {t), (ii) 

which, upon taking the Laplace transform j(p) = 
C[J(t)], g(p) = C[G(t)], gives the expression for the 
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source term, 



hip) 



4A(p-p 2 ) 



Kl+4A 5(0il) -(l-2r)Aft)' 



where we have defined h(p) = 5(o,o)(p) + 5(0,2) (p) + 
(p) f° r the square lattice. 

Note that in the classical voter model (r — 1/2), 
the sink terms vanish (j e — 0), and a single source is 
located at the origin. The calculation may then pro- 
ceed exactly as described in Ref. [22]. On the other 
hand, for the infinite-temperature Kawasaki dynamics 
(r = 0), the source and sink terms create no net correla- 
tion (j e = — jo/4), and they serve only to maintain the 
stationary state Co = p, C-^q = p 2 . 

The discussion so far has been exact. We now turn 
to the long-time asymptotic solution of the correlation 
function, for which it is sufficient to consider the be- 
haviour of jo, j e at small p (p <C A). In the following we 
shall make use of the following expansions: 



a. The density of interfaces cab between A and B 
cells drops asymptotically as cab ~ l/2rlni. 
(12) The density of interfaces between adjacent A and B cells, 
cab(£) = 2(p — C e (t)), is an order parameter used to 
describe the transition from an uncorrelated initial state 
with cab = 2(p — p 2 ), to the jammed absorbing state 
cab(^ — * °°) — 0. For r — 1/2, it was previously shown 
that cab ~ 1/ In t. 

From its definition, we identify the order parameter to 
be proportional to the source term, viz. cab(^) = Jo/2A, 
so that the long-time asymptotic behaviour is found from 
Eq. ( Jl4[ ). This expression reveals the continuous transi- 
tion between Kawasaki dynamics (r = 0, cab = const) 
and voter dynamics (r — 1/2, cab ~ 1/lnt). In par- 
ticular, as t — > oo, the absorbing-state phase transition 
occurs at r = 0. 



lim O(o,o) (p) 

p/A^0 



4ttA 



ln(32A/p) + C[pln(p)], 



and (for the same limit p/X — > 0) 

5(o,o) (p) -5(o,i) (p) = ^ + 0\pln(p)] 
5(0,0) (p) -5(1,1) (p) = ^+0[p\n{p)] 
5(0,0) (p) -5(0,2) (p) = ^ + 0[plu(p)} . 

With these expansions, the source terms take the long- 
time asymptotic values 



P /™o J0(P> p[7r(l-2r)+2rln(32A/p)] 
The corresponding long-time behaviour is therefore 



lim J (t) = — — 

t>l/A 7T(1 



2r) +2rln(32At) 



(14) 



giving the expected 1/lnt decay. Now the effect of 
Kawasaki dynamics becomes clear: one may see that 
reducing r from its maximal value of 1/2 has the op- 
posing effect of weakening the magnitude of the net 
source (Jo — 4J e oc 2r) while extending the time over 
which the source decays, t r ~ (t r =i/2) 1 ^ 2r ■ In the trivial 
limit r = 0, the system becomes stationary as expected, 
whereas when r = 1/2 one retrieves the same expression 
as in Ref. [2"2] . 

Taken together, Eqs. d9|, (fTTl and (Tl4j) give the 



solution for the long-time asymptotic behaviour of the 
two-point correlation function. With a mind towards 
experiment, as well as to compare with the known 
results for r = 1/2, we now summarise the features of 
this solution: 



b. The spatial correlation function Ci(t) decays as 
a(t) — b(t)]n\i\ at short distances, and as a Gaussian 
at long distances. 

Away from the origin, where a continuum description 
suffices, then the correlation function depends only on 
the distance x = and all sources appear to be located 
at the origin, viz. j/ off '- ) = <5 i0 ( Jo + 4J e ). Replacing 
the discrete problem with a continuous one simplifies 
the Green's function, with lim^! Gi(t) — > G(x,t) = 
e -x /4At^(4 7r At) for a square lattice. As a consequence, 
the long-time asymptotic correlation function may be 
approximated as lini t> i/A:r»i Ci(i) = C{x,t), with 



C{x,t)=p 2 



2r{p-p 2 ) 



7r(l - 2r) + 2rln(32At)' 



-Ei 



x^\ 

Ixt) 



( 15 ) 

(13) where Ei(ar) = — j_ x dt e 1 jt is the exponential integral. 
From the small-argument expansion of this integral, we 
find 



lim 

t»l/A, x 2 <At 



C(x,t) 



p- p 



a(t) — b(t) Inx . 



with 



a(t) = 1 - 



1 -2r(7r-7 e -ln8) 
2rln(32At) + 7r(l-2r) ' 



b(t) = 4r/[2rln(32At) +tt(1- 2r)], and j e is the Euler 
constant. 

At long times, C(l,f) > C e (t), so one may infer that 
the correlation function is concave near the origin. Away 
from the origin, it is useful to define a correlation length 
£ to characterise the short-range correlations, viz. 



r^iP-P 2 ) 



dC 
Ox 



= b, 



(16) 



corresponding to the typical size of A-cell rich domains 
growing as £ ~ In t. One may see that when r > 0, a 
variation in r merely adjusts the correlation length by a 
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constant, in contrast to its effect on the order parameter 
cab- For r = 0, the lattice configuration is random and 



C(x,t) 
length. 



= p , corresponding to an "infinite" correlation 



c. The time to saturation T of a lattice of N sites is 
approximately T ~ 7V(lniV+ l/(2r))/A. 
As mentioned earlier, the classical voter model predicts 
that any finite-sized system inevitably approaches an ab- 
sorbing state, in which the lattice is either completely 
saturated by type A cells or else they have become ex- 
tinct. Fortunately, the time-scale T in which the ab- 
sorbing state is reached is T ~ (Nh\N)/X for r = 1/2, 
which, for a biological system with 1/A ~ 1 week and 
N 3> 1, far exceeds the lifetime of a mammalian organ- 
ism. 

For arbitrary r, the time-scale T may be estimated 
by repeating the calculation in Ref. [TH] : The saturation 
condition is Ci(T) = N p. Replacing the summation 
by integration, Ci — > xC(x, t)dx, we arrive at the 
asymptotic relation 



lim > ' 



d- p 2 ^ ArXT 
p- p 2 ~ tt(1 -2r) + 2rln(32AT) ' 



and the result for T follows. Interestingly, this result 
predicts that in even a reasonably large system and fi- 
nite r, when N 3> e n l 2r , then the time to saturation is 
insensitive to the value of r. Indeed, for such systems 
the cross-over to non- voter-like behaviour only occurs at 
very small values of r ~ 1/ In TV. 



d. The product of the correlation length and interface 
density gives a time-invariant characteristic of voter-like 
coarsening. 

It has been previously noted that in voter-like coarsen- 
ing, the characteristic length-scale of domains is inversely 
proportional to the interface density cab |16j . An im- 
portant implication of this observation is that at asymp- 
totically long times, one may identify a time-invariant 
characteristic of the correlations, which we define as 
n = (p- p 2 )/(c A bO = 2r/7r. 

To characterise the fi-constant, one may identify its 
definition as the ratio between the domain size £ and 
domain circumference, as calculated from the number 
of A-B interfaces associated with each domain cab£ 2 - 
Thus, the domain perimeter has a trivial fractal dimen- 
sion of one, with f2 indicating the perimeter roughness, 
or curvature. Values of f2 ~ 1 may be associated with 
cohesive and smooth domains (i.e. with a large area- 
to-interface ratio), whereas systems with O — ^ have 
highly fragmented, or rough, domains. Not surprisingly, 
smaller values of r lead to rougher domains as a result of 
the Kawasaki dynamics. Yet, one may see that even at 
the maximum value of r — 1/2, the voter model predicts 
rough domains (Q < 1/2) — an observation readily seen 
in simulations, see e.g. Ref. |23j . 



With an eye to the empirical analysis in the next 
section, let us note that the fi-constant allows us to 
characterise static observations of type A cell correla- 
tions, and thus provides a valuable test of whether a 
given data set is consistent with the long-time behaviour 
of process Q . 

e. The lattice geometry influences only the timescale 
of clustering, with higher coordination number corre- 
sponding to a faster timescale. 

For lattice geometries with different coordination num- 
ber z 7^ 4 (and one site per unit cell), the calculation 
proceeds as for a square lattice, but now three modifi- 
cations must be made: First, as mentioned above, the 
central source term is J = zX{p — C e ), and there are 
now z sink terms J e at nearest neighbouring sites. Sec- 
ond, the appropriate Green's function now has the (non- 
separable) form 



G,(t) 



1 

4^2 



d qexp 



zq • l 



Xt 



«q e 



from which we obtain the general form of the small- 
p (long-time) expansion of the Laplace transform 
lim p / A ^ g Q (p) ~ ln(X/p) + d /2X and ,g ; (p) = g (p) - 
d[/2X, where d\ are numerical constants, e.g. d e = 2/z. 



Third, the sum over sink terms in Eq. (10 1 is revised to 
reflect the lattice geometry. With these three modifica- 
tions, one finds that the order parameter cab takes the 
general form cab = 27r(p — p 2 )/[a(l — 2r) + 2rln(/?Ai)], 
where a and (3 are geometry-dependent constants. One 
is led to conclude that the lattice geometry only serves 
to rescale the time variable by some constant, t — ► 
[3e atyl ~ 2r H. It is interesting to contrast this with the 
more significant effect of modifying the branching ratio 
r, which instead rescales hit. 

With regards to the experimental system, the appar- 
ent insensitivity of the results to details of the lattice 
geometry reinforces the validity of a lattice-based de- 
scription of the basal layer. Namely, while it is clear 
that the basal-layer is not a periodic lattice of uniform 
cell size and coordination number, it may nonetheless be 
modeled as such. 

Finally, because the average coordination number of 
z = 6 is expected for the biological system, we have 
calculated the order parameter exactly for a hexagonal 



lattice, giving a 



(hex) 



0.98tt, and /3< hcx ) = 51. 



This completes our theoretical discussion of process Q , 
and of the monomer-monomer surface catalysis model. 
It now remains to be seen whether the exact model of 
cell division Q is indeed described by the properties of 
the approximate process Q (with uniform A). Let us 
recall that the latter is a reasonable approximation as- 
suming that the density of type A cells is approximately 
uniform. Thus, although the approximation must fail on 
the time-scale XT ~ N In N associated with the jamming 
transition, the logarithmically slow growth of correla- 
tions suggests that at shorter times (t <C T) the degree 
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FIG. 3: Comparison of the exact basal layer lattice model |2| 
(red curves) with the simplified model Q (black), (a) The 
order parameter 1/cab is plotted against Int. The curves cor- 
respond to the simplified model with r = 1/2 (dashed black) 
and r — 1/4 (dotted), and to the exact model with r — 1/2 
(dashed red) and r — 0.08 (solid curve). For the simplified 
model, the analytical expression for cab was used. For the 
exact model, the results were obtained from numerical sim- 
ulation on hexagonal lattices of size N = 1024 x 1024, with 
p = 0.22 and setting A = 1. Inset: The correlation function 
plotted against logarithm of the distance In a;, evaluated from 
numerical simulation and Eq. (171. The time-invariant ratio 
(p - C(x, *))/cab is plotted at At = 50, 100, 1000 for the clas- 
sical voter model (dashed), and at At — 350,3500,7000 for 
the exact model with r = 1/2 (solid curves). Consistent with 



the expected behaviour (see Eq. ( 15 1), the curves overlap and 
are convex near x — 1, becoming linear at x > 1. (b) The 
long-time asymptotic roughness constant SI plotted against r 
for both models. Data points correspond to the values calcu- 
lated using Eq. [17] from numerical simulations such as shown 
in Fig. [4] using the algorithm described in section |III A| Er- 
ror bars result from fluctuations due to finite-size effects, as 
evaluated by considering the random variation in f2 over the 
time of the simulation. The black curve gives the theoretical 
value of Q, for the simplified (monomer) model on a hexagonal 
lattice, while the red curve gives the best-fit to the numerical 
results for the exact model. 



of progenitor cell clustering should be sufficiently low as 
to make the analysis self-consistent. 



To test whether the model indeed satisfies the ex- 
pected behaviour, in Fig. [3] we compare the spatial cor- 
relation of type A cells, as given by the properties de- 
rived above, with the results of process ^ as obtained 
by numerical simulation using a Gillespie-like algorithm 
(described below in section III A). An example of the 
cellular automata simulations used for the comparison is 
shown in Fig. [4] where the distribution of type A cells 
(black) is shown on a hexagonal lattice. Qualitatively, 
one may see from this figure that process Q does indeed 
give rise to some clustering. By extracting the correla- 
tion function from such figures, we obtained the quanti- 
tative comparison shown in Fig. [3] In Fig. [3^, we com- 
pare the evolution of the inverse order parameter 1 / cab 
against lnt for the two models. One may see from the 
linear behaviour of the exact model (red online), that 
the order parameter indeed shows the expected 1/ In t 
decay. Equally, we may confirm that the radial correla- 
tion function C(x,t) has the functional form a — b In x, 
with parameters (1 — a), b both proportiaonal to cab> 
as demonstrated by plotting (C(x, t) — p) /cab in Fig. [3^ 
(inset). Finally, we may confirm that both models have 
the same long-time r-dependence by plotting f2 against 
r in Fig. [3J). Here, one may see that indeed fl depends 
linearly on r. 



Yet, although Fig. [3] reveals that processes ^ and 
Q result in the same functional dependence of the cor- 
relation function on r and t, it is striking that the 
exact model ^ results in a slower growth in correla- 
tions compared to the simplified model, characterised 
both by a slower decay of cab seen in Fig. [3k and by 
rougher domains (or lower fi) seen in Fig. |3b. How 
might one explain this change? Referring to the dis- 



cussion in section II A let us recall that the simpli- 
fied model is connected to the exact model by relating 
the local division rate to the mean-field vacancy density 
(A(x) = A'n (x)), with the two models becoming equiv- 
alent when 7i0 (x) = const for an arbitrarily small degree 
of coarse-graining. This condition is satisfied in the limit 
tr > A' with the process A — > A removed. How- 
ever, for the physically meaningful case a <C A' one can 
no longer treat the vacancy density as uniform, leading 
to quantitative (but not qualitative) differences in the 
behaviour of the two models. As shown schematically 
in Fig. [5j the vacancy population has the effect of ac- 
celerating cell division on the edge of smooth progenitor 
cell clusters, as a result of the higher concentration of 
nearby post-mitotic cells migrating into the super-basal 
layers. Conversely, cell division on rough cluster edges is 
slowed down. In total, rough A-B interfaces remain sta- 
ble over a longer period of time, leading to the observed 
differences between the two models. 
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FIG. 4: Cellular automaton simulation of process (|2|, show- 
ing the distribution of progenitor cells (black hexagons) on a 
lattice of N = 200 x 200, and using the experimental branch- 
ing ratio r = 0.08 and progenitor cell fraction p = 0.22. The 
frame shown corresponds to an evolution time of t = 30/A, 
where t = corresponds to random initial conditions. White 
areas are fully occupied by post-mitotic cells. 




FIG. 5: (Color online) Schematic demonstrating the varia- 
tion in the effective division rate in process accounting 
for the quantitative variation in behaviour between processes 

Sand Q. The reactive interface between progenitor cells 
;ht grey, red online) and post-mitotic cells (dark grey, blue 
online) is indicated by the thick grey line, with a faster divi- 
sion rate shown in light grey and slower division rate shown 
in dark grey. On the left, a smooth interface results in faster 
cell division, as few progenitors on the boundary must com- 
pensate for the migration of many post-mitotic cells in the 
bulk. On the right, the "rough" interface results in a vari- 
ation between fast and slow division rates, due to restricted 
access of vacancies into the rough interfacial regions. Refer- 
ring to the discussion at the end of section [il B[ this leads to 
increased stability of rough compared to smooth interfaces. 



III. EMPIRICAL ANALYSIS 

We are now in a position to turn to the experimental 
analysis of the IFE. In order to identify the underlying 
rules of cell division and differentiation, one can envision 
two types of experiment: 

• First, one may track the fate of individual cells 
and their progeny, and then look for a cell ki- 
netic description compatible with their observed 
behaviour. Such an approach was used with con- 
siderable success in the definition of the zero- 
dimensional process ([I]) through clonal analysis, 
and we shall extend it to the analysis of the spatial 
process in section |III A| 

• Second, in section [TlIB| we shall consider the statis- 
tics of the entire population of basal layer cells. In 
particular, by immunostaining basal layer cells for 
markers of cell proliferation, it is possible to anal- 
yse the spatial distribution of all progenitor cells 
within the layer. Then, referring to the theoretical 
discussion in section [TTJ one may look for a signa- 
ture of the underlying cell kinetics in the spatial 
correlation of proliferating cells. 

Note that the two types of experiment give access 
to independent aspects of cell behaviour: The former 
probes the temporal evolution of cell lineages, whereas 
the latter reveals the static basal layer morphology. As 
such, the experiments provide a significant degree of mu- 
tual verification of any proposed theory of cell behaviour. 
Yet, even the best of such experiments leave room for 
some ambiguity: For example, it is far from clear what 
importance should be assigned to the embryonic devel- 
opment of the IFE in predetermining the spatial distri- 
bution of cells. Moreover, one may in principle conceive 
of regulatory pathways that leave no signature on the 
spatial distribution of cells, or which may not be distin- 
guished from an independent stochastic process. Thus, 
in the following we will look for the simplest possible 
model of cell behaviour that succeeds in capturing the 
known biological constraints. 

A. Clonal analysis 

To begin our analysis of the empirical data, we start 
by considering the fate of individual labelled cells and 
their progeny, hereafter referred to as clonal fate data. 
In an extensive experiment reported elsewhere [5], the 
low-frequency labelling of approximately 1 in 600 basal- 
layer epidermal cells at a defined time was achieved by 
a drug-induciblc genetic event, which resulted in expres- 
sion of the enhanced Yellow Fluorescent Protein (EYFP) 
gene in a cohort of mice. At intervals, the EYFP label 
was detected by confocal microscopy, which enables 3D 
imaging of entire sheets of epidermis. By analysing sam- 
ples of mouse epidermis at different time points, it was 
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FIG. 6: (Color online) (a) Comparison of the empirical clone 
size distribution (data points) to predictions of process |2| 
(solid curves), as obtained from Monte-Carlo simulations of 
10 4 labelled clones. The distributions are plotted in terms of 
the probability Vk{t) for a clone to have between 2 fe_1 + 1 
and 2 k basal layer cells at time t after labelling, normalised 
to include only clones with 2 or more cells in the basal layer 
(k > 1). The empirical data is reproduced from Ref. [5]- (b) 
Examples of the basal layer structure of several large late- 
stage clones evolving according to process (pi), starting from 
a uniform random distribution of unlabelled cells, with one 
single type-A cell labelled at t = 0. Light and dark grey 
hexagons (red and blue online) indicate sites occupied by 
type A and B cells, respectively. White areas are populated 
by unlabelled cells. Each frame corresponds to the progeny 
of one initially labelled cell. 



possible to analyse the fate of labelled basal layer clones 
at single cell resolution in vivo for times up to one year 
post-labelling in the epidermis (see, for example, Fig.JTJ)) 

mm- 

The distribution of cells within a clone constitutes a 
valuable data set, which one may analyse for signatures 
of the underlying rules of cell division. By scoring the 
total number of labelled cells n in each clone at progres- 
sive time points, it was possible to access the evolution 
of the full clone size distribution, as given by the proba- 



bility P n >o(t) f° r finding a clone with n > 1 basal layer 
cells at a time t post-labelling, shown in Fig. ga). The 
properties of the size distribution were used to infer the 
laws summarised in process 

@ [5, 6J. In particular, it 
was shown that in the long-time limit (t ^S> \/r\), the 
basal layer clone size distribution conforms to the scaling 
form 



lim P n>0 (t) 



r\t 



exp 



\rXtJ 



As discussed in section |TTJ we must first establish, for 
any proposed model of spatial behaviour, that the zero- 
dimensional clone size distributions are faithfully repro- 
duced. 

To this end, we conducted multiple simulations of pro- 
cess ([2]) as an asynchronous cellular automata evolving 
on a hexagonal lattice of N — 60 x 60 sites over a pe- 
riod corresponding to T = 60 weeks in the experimental 
system (recall that A = 1.1 /week [5]). At t = 0, the lat- 
tice was fully occupied by randomly-placed type A and 
type B cells. A single, randomly chosen, type A-cell was 
assigned a hereditary 'label' at the start of each simula- 
tion. In effect, each such simulation mimics the evolution 
of one labelled clone, so that repeated simulations may 
be used to sample the full clonal statistics, viz. Montc- 
Caiio sampling. In particular, by tracking the number 
of labelled cells as a function of time over 10 4 such clone 
simulations, we could compare the clonal statistics pre- 
dicted by process ([2| with those expected from the zero- 
dimensional process ([T]). In keeping with the empirical 
fit from the previous section, we used a B-ccll migration 
rate of T = 0.31/week, and an initial A-cell fraction of 
p = 0.22. The 'fast' rates of cell division and hole dif- 
fusion were set to be A' = 10 4 /week and a — 200/week, 
although the precise values are unimportant provided 
that A' > a > T, as discussed in section [Tl] The results 
are plotted in Fig. pla) , where we compare the clone size 
distribution P n >o(t) to the empirical data. One may see 
that the fit is remarkably good, and, within the current 
empirical resolution is indistinguishable from the predic- 
tions of zero-dimensional model (cf. Ref. |6J. In sum- 
mary, the fit provides a first validation of process ^ as 
a viable model of spatial behaviour in the basal layer. 

We may now revisit the clone fate data with a view 
to study spatial structure. Unfortunately, although the 
clone size data was stored for all clones, images show- 
ing their spatial structure were retained only for a small 
number of clone samples. Therefore we are not in a 
position to conduct a comprehensive quantitative anal- 
ysis of clone shape evolution. Nevertheless, a striking 
qualitative feature of clones is that they remain largely 
cohesive, as shown in the example in Fig. [T 3. (Indeed, 
without this property the very enterprise of clonal analy- 
sis would have proved difficult). Therefore, to challenge 
the validity of the proposed lattice model of cell divi- 
sion, we test, using additional Monte-Carlo simulations, 
the ability of the model to produce cohesive clones over 
the one- year time period of the experiment. Represen- 
tative results are shown in Fig. [6] for large clones at 60 
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weeks post-labelling. When A' ^> a ^> T, clones were 
seen to remain largely cohesive throughout the period 
of the simulation. We may therefore conclude that, at 
least with respect to the existing clonal fate data, pro- 
cess Q presents a reasonable phenomenological descrip- 
tion of the spatial behaviour of basal layer cells. In par- 
ticular, the cohesive nature of clones may be completely 
explained in terms of an independent stochastic process, 
with no evidence for further forms of regulation. 

Yet, to what extent is the observation of cohesive- 
ness sensitive to modifications of the basal layer lattice 
model? To test the degree to which cohcsivcncss con- 
strains the model, we allowed for a degree of cell mo- 
bility within the basal layer, by introducing the addi- 
tional exchange process AB — > BA with rate 7 ~ V. The 
exchange process is a reasonable candidate for cell be- 
haviour, motivated by the observation that keratinocytes 
in culture are highly motile, which leads one to postulate 
whether cells in vivo are capable of independent lateral 
migration in the basal layer. We found that for any non- 
small value of 7, (i.e. 7 > T), the progeny of labelled 
clones rapidly dispersed, as shown in Fig. [7] We are led 
to conclude that epidermal cells in vivo move only in re- 
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FIG. 7: (Color online) The loss of cohesiveness with increas- 
ing the relative hopping rate a/X' is shown through exam- 
ples of late-stage clone simulations. Light and dark grey 
hexagons (red and blue online) indicate sites occupied by type 
A and B cells, respectively, from the same clone. White ar- 
eas are populated by unlabelled cells. The simulations used 
the same parameter set as described above for the clones 
in Fig. pi but using the parameter values: a = 20/week, 
a = 200/week with A' = 10, 000/week for the top-left and 
top-right clones respectively; a = 200/week, a = 2000/week 
with A' = 200/week for the bottom-left and bottom-right 
clones. For the latter, which is clearly unphysical, the scale 
has been reduced twofold to demonstrate the wide dispersion 
of labelled cells. The examples demonstrate that clones are 
cohesive when a/X' <C 1, but dispersive otherwise. 



sponse to a local density gradient resulting from cell divi- 
sion and migration. Moreover, a second investigation in 
which the hard-core mobility a was made larger than A' 
(and keeping 7 = 0) again led to a loss of clone cohesive- 
ness. We are therefore led to conclude that, within the 
framework of the non-interacting lattice model, progen- 
itor cells undergo division rather than lateral migration 
as a response to a local drop in density. 

In conclusion, the observation of clone cohesiveness 
imposes a severe constraint on cell behaviour, which al- 
lows us to rule out several variations of the basic lattice 
model. Yet, at least qualitatively, there is no clear evi- 
dence either for or against additional regulation of cell 
division in the clone shape data. So to what extent is 
the lattice model truly capable of shedding new light on 
the mechanism of cell fate regulation? To do better, we 
must look for quantitative data for comparison. Fortu- 
nately, such data can be found in the readily-available 
spatial distribution of progenitor cells. 



B. Correlation analysis 

As mentioned earlier, an intriguing feature of the non- 
interacting lattice model is the prediction of clustering 
of proliferating cells. How does this prediction compare 
with experiment? Before we turn to consider new re- 
sults for mouse tail-skin, it is interesting to first consider 
results presented in past work. In a detailed study of pro- 
liferating cells in Hamster cheek epidermis by Gibbs and 
Casarett j9], a subset of progenitor cells were labelled 
using a radioactive marker for DNA synthesis (S-phase). 
By measuring the number of unlabeled cells separating 
consecutive labeled ones in one-dimensional basal layer 
cross-sections, it was possible to access the full radial 
distribution for the separation between adjacent cells un- 
dergoing S-phase. Remarkably, while the large-interval 
distribution decayed exponentially as expected for an un- 
correlated random distribution, the probability of find- 
ing a nearby cell in S-phase was significantly higher at 
short distances, indicating that proliferating cells were 
clustered. A second study addressing the distribution of 
S-phase cells in mouse esophagus also revealed identical 
qualitative results [TO] . 

Not surprisingly, in the absence of the intuition af- 
forded by the voter model, such clustering has been in- 
terpreted in the biological community as evidence of an 
underlying regulatory process, which leads to the syn- 
chronous division of nearby cells [10 . On the other hand, 
an independent analysis of clone fate data in Ref. [5J in- 
dicates that cell division within clones occurs indepen- 
dently It is therefore satisfying to note that the ten- 
dency of proliferating cells to cluster is in fact consistent 
with independent division — indeed it is the hallmark of 
voter-model dynamics. 

To extract the spatial correlation between progeni- 
tor cells in mouse tail skin, we analysed confocal mi- 
crographs of basal layer cross-sections of IFE that were 
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FIG. 8: (Color online) (a) Confocal micrograph of whole- 
mounted mouse tail skin IFE, showing the two-dimensional 
basal layer immunostained for the nuclear marker DAPI 
(dark grey, or blue online) and the proliferation marker Ki67 
(light grey, or red online). The area surrounding the stained 
nuclei is occupied by unstained cell cytoplasm (black), (b) 
The empirical radial correlation function C(x), as defined in 
Eq. (18 1. Data points show results obtained by analysis of 



the Ki67-stained epidermal wholemounts exemplified in (a), 
taken from a mouse aged 8 weeks; the dashed line shows the 
fit to the analytical form of the correlation function p red i cted 
by process m, C(x,t) = a(t) -£~ 1 (t)\nx (see Eqs. 15 16 1, 
where the single time point t is fixed by the age of the mouse 
and by initial conditions (see discussion in section [rV| . From 
the fitted slope one may extract the correlation length £(£) 
defined in Eq. (16 1, giving £ = 14.7 cell diameters. 



stained for the proliferation marker Ki67, such as shown 
in Fig. |8ja). Cells bright in Ki67 are designated as type 
A (progenitor) cells, whereas Ki67-dull cells were desig- 
nated as type B (i.e. differentiated) cells, with no capac- 
ity to divide. Using image analysis software (Image J), 
the coordinates of each Ki67-bright cell were extracted. 
This data allows a full statistical analysis of the spatial 
distribution of progenitor cells. In particular, we shall 



focus on the radial correlation function 

^ de 



C(x) 



27T J A 



dx'n(x')n(x' + x(a;,0)) ) (17) 



where A denotes the area of each sample, n(x) denotes 
the areal density of proliferating cells at position x, and 
the brackets (•) indicate averaging over all basal layer 
samples. 



Using Eq. (17 1, the aim of the correlation analysis is 



to assess whether basal layer progenitor cells in adult 
mice do indeed cluster, and if so, to assess whether the 
experimental correlation function is consistent with the 



predictions made in section II B In principle, one may 
look to the data for signatures of the expected spatial 
dependence C(x,t) ~ a(t) — b(t)\nx (for a given value 
of t) , as well as for evidence of increased clustering with 
time, viz. cab(^) ~ 1/(27- In i). For the latter, however, 
the temporal analysis is difficult to implement due to 
the sensitivity required to resolve the 1/ In t decrease in 
cab at long times. In particular, the predicted increase 
in clustering over a biologically-relevant period of Xt ~ 
10 1 — 10 2 cell cycles corresponds to a decrease in cab of < 
10%, whereas variations in the efficiency of Ki67 labelling 
between different mice introduce systematic errors of the 
same order. Therefore, in the following we shall restrict 
ourselves to the quantitative analysis of tissue samples 
taken from a single adult mouse (aged 8 weeks). 

To incorporate the empirical coordinates 
of proliferating cells into Eq. (17), 
place the product n(x')n(x' + x) with 
J2i j <5( x ' — Rj)/(x' + x — Rj), where 6 (it) is the Dirac 
delta function, {R^} is the set of all progenitor cell 
coordinates, and /(X) is a two-dimensional Gaussian 
envelope with width w. With this substitution, the 
integrals in Eq. (17 1 can be solved exactly. Then, 



we re- 
the sum 



explicitly accounting for the averaging procedure over 
sample of variable size, we obtain the expression 



C(x) 



2-kw 2 N a (x) 



N A N A (x) / 2 , p2 " 

E E «p 3 



2w 2 



xln 



xRi 



(18) 



Here, the sum i, j over all progenitor cell coordinates 
arises from the empirical expression for n(x')n(x' + x) 
given above, and we have defined Rij = |Rj — Rj|. The 
prefactor, exponential factor and modified Bcssel func- 



tion (Iq) result from solving the integrals in Eq. (17 1 



The total number of progenitor cells found in each sam- 
ple is Na, and Na(x) is the number of progenitor cells 
at a distance x or more from the sample edges. As 
defined, the correlation function avoids errors result- 
ing from edge effects by averaging each sample over 
the Na(x) progenitor cells that are unaffected by the 
finite sample size. To correctly average over differ- 
ent samples, as indicated by (•), the sample results are 
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weighted by N A (x), e.g. C {1+2 \x) = [N^ ] {x)C {1 \x) + 
Nf{x)C^{x)]/[N^{x)_+Nf{x)]. 



Making use of Eq. (18 1, the experimental correlation 



between progenitor cells was averaged over 10 samples of 
approximately 30 x 30 cells each, see Fig.|8jb). Remark- 
ably, the data indeed shows a significant degree of pro- 
genitor cell clustering, in good agreement with the lin- 
ear decay in In a; expected from Eq. (151. In real terms, 



one may infer from the value of the nearest-neighbour 
correlation C(l) that each progenitor cell is in contact, 
on average, with approximately two adjacent progenitor 
cells, compared to 1.3 progenitor cells expected for an 
uncorrelated random distribution (at p = 0.22). 

For a more careful test of the theory, one may extract 
from Fig. [8jb) the order parameter cab = 0.29 ± 0.02, 
and the correlation length of £ = 14.7 ± 0.7 cells, 
from which we find the empirical roughness constant, 
fi = 0.04 ± 0.01. Repeating the analysis with samples 
taken from different mice results in the same value of 
cab, but with values of f2 in the range = 0.02 — 0.04. 
Referring to Fig. |3jb), where the empirical value of Q 
is compared to the model predictions (dashed), we see 
that the model indeed recovers the correct order of mag- 
nitude of the roughness constant. Qualitatively, then, it 
appears that the model is consistent with the observed 
clustering. Let us emphasise that this fit requires no ad- 
ditional parameters, and is purely a result of mapping 
the cell kinetics onto a lattice. 



IV. DISCUSSION AND CONCLUSIONS 

To summarise, we have demonstrated that the zero- 
dimensional model of cell division is consistent with the 
maintenance of the basal layer at uniform density, and 
we have shown that the size distribution and qualitative 
shapes of labelled clones are consistent with a simple 
stochastic model of cell division and differentiation on a 
two-dimensional lattice. These results explain why the 
recently-discovered "zero-dimensional" model of cell be- 
haviour gives such an excellent fit to the clone fate data 
despite taking no account of additional regulatory path- 
ways. Significantly, despite the many forms of cell fate 
regulation known to exist in development and adult tis- 
sue, the only extra-cellular regulation required to under- 
stand the existing observations of clone fate in normal 
IFE is steric, viz. the coupling of cell division to the 
local cell density. 

Beyond the success of the model in explaining the 
observed clone fate data, we have also identified that 
the degree of progenitor cell clustering, as measured by 
Ki67 staining, is in good qualitative agreement with the 
predictions of the spatial process. In particular, there 
are two features that allow us to characterise the spa- 
tial process. First, the empirical roughness constant 
O = 0.04 ± 0.01 has the expected order of magnitude 
predicted by the model (Fig. |3]d), and second, the cor- 
relation function is in excellent agreement with the ex- 



pected decay form a(t) — b(t) \nx (at fixed t), as seen in 
Fig-i 

Taken together, these results have two important 
implications for future investigations of epidermal cell 
fate regulation: First, by demonstrating that the zero- 
dimensional process ([!]) is indeed capable of maintaining 
a uniform total basal layer cell density, the spatial pro- 
cess consolidates the proposed stochastic model as a ro- 
bust platform for investigating biochemical constituents 
in future work. For example, by over/under-expressing 
specific genes and then studying the resulting change in 
the empirical parameters (r, A, p), one may attempt 
to identify the role of each constituent in regulating 
cell behaviour. Second, the new model introduces a set 
of spatial measures (such as the roughness constant fi) 
that yield further information in such investigations of 
the biochemistry, beyond that which may be obtained 
through empirical evaluation of the zero-dimensional pa- 
rameters alone. 

Beyond the qualitative features of the dynamics, one 
may ask whether there are any implications to the quan- 
titative features of the correlation analysis. Unfortu- 
nately, as mentioned earlier, variations in the efficiency 
of Ki67 labelling prevent us from generating the compre- 
hensive statistics necessary to accurately quantify the 
universal values of the system parameters. Neverthe- 
less, taken at face value, it appears that the empiri- 
cal value of £1 obtained from the sample analysed in 
section IIIB is only consistent with process ^ when 
one imposes a branching ratio r — 0.19 ± 0.04, which 
differs significantly from the value of r = 0.08 estab- 
lished through clonal analysis (see fig. [3b) . This value is 



further consistent with the observed AB interface con- 
centration of cab = 0.29 and the correlation length of 
£ — 14.7 cells, which correspond to a long evolution time 
of Xt w 10 4 — 10 6 for r = 0.08, but a realistic biological 
time-scale of \t « 10 1 - 10 2 for r = 0.19. 

Although it is possible that the discrepancy in the in- 
ferred value of r results from errors in progenitor cell 
classification as discussed in section IIIIBI the difference 



is large enough to call into question the reliability of the 
correlation analysis. Beyond the issue of progenitor cell 
labelling efficiency, there is also the more basic question 
of whether Ki67 is at all effective as a marker of progen- 
itor cells. In particular, it is widely accepted that Ki67 
is a marker of cell growth as well as proliferation, which 
may imply that differentiated cells remain Ki67-bright 
for some time after division, or that progenitor cells may 
fail to continuously express Ki67 [24j|25]. However, as- 
suming that Ki67 is indeed a faulty marker, it becomes 
difficult to explain why the analysis nevertheless results 
in an excellent fit to the 'a(i) — b(t) In x' decay form of 
correlations (at fixed t). 

We may also challenge our assumption of using "ran- 
dom initial conditions" to model steady-state mainte- 
nance. If initial conditions are at all relevant, then the 
observed clustering may be a signature of development 
and growth processes that are no longer active in the 
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adult system. To test the importance of initial condi- 
tions, we have compared the order parameter cab be- 
tween adult mice aged 8 and 60 weeks (corresponding 
to yound and old adult mice). If clustering is inherent 
to initial conditions, then one would expect cab to grow 
over time as the memory of initial conditions erodes. On 
the other hand, as mentioned in section |IIIB[ steady- 
state maintenance results in < 10% change in the order 
parameter over this time period, which would be unde- 
tectable given the systematic errors in Ki67 labelling. 
Indeed, we find no change from cab = 0.29 ± 0.02 at 60 
weeks, so there is no indication that the observed clus- 
tering is a signature of tissue development. 

Finally, it is interesting to speculate whether addi- 
tional forms of regulation may be capable of reconcil- 
ing the higher value of r = 0.19 found for the spatial 
process, with the lower value of r = 0.08 found from 
the zero-dimensional analysis. We propose a simple re- 



vision of the model that is capable of reconciling the 
spatial and zero-dimensional data: Referring to process 
2 ) , if the two channels of asymmetric division in process 
2 ) occur with different probabilities (viz. Pa0^ab =/= 
-Pa0-+ba ^ 1/2 - r, and Pa0^ab + Pa^ba = 1 - 2r), 
then the effective value of the branching ratio r is effec- 
tively renormalised in the spatial process, while leaving 
it unchanged in zero-dimensions. To reconcile the two 
empirical values of r, it is simple to show that one re- 
quires Pao^ab — 7Pa0^ba — (1 - 2r)/8. Thus, one 
might speculate that the increase in clustering is asso- 
ciated with a spatial asymmetry in daughter cell fate 
during asymmetric division. It would be an interesting 
challenge to devise an experiment with which to test this 
asymmetry. 
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